;;....................................................................................
;; Updraft velocity PDF
;; Kai Zhang, PNNL 
;;....................................................................................
   load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
   load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
   load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
   load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"   


   yes_incloud = True
   yes_all     = False

   yes_incloud = False  
   yes_all     = True

   fna = "wsub3h_sp_T63L31K01_2010_jan-jun.nc" 
   fnb = "icnc3h_sp_T63L31K01_2000_jan-jun.nc"

   fla = addfile(fna,"r") 
   flb = addfile(fnb,"r") 

   vw = fla->W_TURB
   nw = flb->ICNC
   vx = fla->W_TURB

   if(yes_incloud) then 
      ;;vw = where(vw.gt.0.001.and.nw.gt.1.e3,vw,-999) 
      ;;vx = where(vx.gt.0.001.and.nw.gt.1.e3,vx,-999) 
      vw = where(nw.gt.1.e2,vw,-999) 
      vx = where(nw.gt.1.e2,vx,-999) 
   else
      if(yes_all) then 
      else
         vw = where(vw.gt.0.001.and.nw.lt.1.,vw,-999) 
         vx = where(vx.gt.0.001.and.nw.lt.1.,vx,-999) 
      end if 
   end if 

   vw = where(vw.le.3.0,vw,-999) 
   vx = where(vx.le.3.0,vx,-999) 

   vw@_FillValue = -999 
   vx@_FillValue = -999 

   vwh = dim_avg_Wrap(dim_avg_Wrap(dim_avg_Wrap(vw(mlev|:,lat|:,lon|:,time|:))))    
   vxh = vx

   pres = 1000.*(/ 5.e-6, 1.e-6, 5.e-5, 1.e-5, 0.000195429078303275, 0.00165527942590415,  \ 
    0.0060569163179025, 0.0147566441446543, 0.0286470074206591,  \ 
    0.0482312496751547, 0.0736913084983826, 0.104949295520782,  \ 
    0.141722559928894, 0.183572381734848, 0.229945927858353,  \ 
    0.280211985111237, 0.333690196275711, 0.38967365026474,  \ 
    0.447445213794708, 0.50628736615181, 0.565485417842865,  \ 
    0.624324411153793, 0.682079493999481, 0.737999826669693,  \ 
    0.791286110877991, 0.841061413288116, 0.886335849761963,  \ 
    0.925964534282684, 0.958599209785461, 0.982633352279663, 0.996140748262405 /) 

             ;;3.64346569404006, 7.59481964632869, 14.3566322512925,                   \       
             ;;24.6122200042009, 38.2682997733355, 54.5954797416925, 72.0124505460262, \  
             ;;87.8212302923203, 103.317126631737, 121.547240763903, 142.994038760662, \  
             ;;168.225079774857, 197.908086702228, 232.828618958592, 273.910816758871, \  
             ;;322.241902351379, 379.100903868675, 445.992574095726, 524.687174707651, \  
             ;;609.778694808483, 691.389430314302, 763.404481112957, 820.858368650079, \  
             ;;859.53476652503, 887.020248919725, 912.644546944648, 936.198398470879,  \ 
             ;;957.485479535535, 976.325407391414, 992.556095123291 /) 
   pres@units = "hPa" 

   hgt = gsn_geop_hgt(pres) 

   print(hgt) 


   hgtx = conform(vxh,hgt,1) 

   print(sprintf("%9.3f",vwh)+"   ") 



   ;va = (/0.49,0.58,0.48,0.62,0.77,0.6,0.425,0.5,0.35/) 
   ;va = va / 1.5 
   va = (/ \ 
               -999, \  
               -999, \ 
          0.4305919, \ 
          0.8094104, \ 
          1.0245100, \  
               -999, \ 
               -999, \
          0.5521003, \ 
          0.5286052, \ 
          0.5390244, \ 
          0.3004771, \ 
          0.5533047, \ 
          0.6982606, \ 
          0.2963851 /) 

   va@_FillValue = -999 

   va = va(::-1) 

  ;ha = (/11.75,11.25,10.75,10.25,9.75,9.25,8.75,8.25,7.75/)  
   ha = (/ 0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,13.5/) 

   ha = ha(::-1) 


   
 plot = new(3,graphic)                          ; create graphical array

;************************************************
; plotting parameters
;************************************************
 wks   = gsn_open_wks ("pdf","xy")                   ; open workstation


 gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")       

 res                   = True                       ; plot mods desired
 res@gsnDraw          = False
 res@gsnFrame         = False

;res@trYReverse        = True                       ; reverse Y-axis
 res@xyDashPatterns    = 0                         ; choose dash patterns

 res@trXMinF           = 0. 
 ;;res@trXMaxF           = 0.7  
 res@trXMaxF           = 3.0
 res@trYMinF           = 7. 
 res@trYMaxF           = 14.

 res@tiXAxisString = "Updraft velocity (m/s)"                   ; x-axis label
 res@tiYAxisString = "Height (km)"              ; y-axis label

 res@gsnStringFontHeightF = 0.025

 res@mpShapeMode = "FreeAspect"
 res@vpWidthF  = 0.6 
 res@vpHeightF = 0.8 
  


 res@xyLineThicknesses = 3.0
 res@xyLineColor     = "black"
 res@gsnCenterString = ""
 ;;plot(0) = gsn_csm_xy (wks,vwh(10:16),hgt(10:16),res) ; create plot
 plot(0) = gsn_csm_xy (wks,vwh(12:18),hgt(12:18),res) ; create plot

 res@xyLineColor     = "red"
 plotx   = gsn_csm_xy (wks,va,ha,res) 

 overlay(plot(0),plotx) 



 opt = True 
 opt@bin_min = 0.001 
 opt@bin_max = 3.0 ;; 0.25  

 res@gsnDraw              = False    ; do not draw picture
 res@gsnFrame             = False    ; do not advance frame
 res@cnFillOn             = True     ; turn on color fill
 res@cnFillMode           = "RasterFill"       ; Raster Mode
 res@cnLinesOn            = False    ; no contour lines
 res@cnLineLabelsOn       = False    ; no contour line labels
 res@cnInfoLabelOn        = False
 res@gsnSpreadColors      = True     ; use full colormap
;res@cnLabelBarEndStyle   = "ExcludeOuterBoxes"
 res@lbLabelAutoStride    = True
 res@tiXAxisOffsetYF      = 0.105

 res@lbOrientation = "vertical"

 res@cnLevelSelectionMode =  "ExplicitLevels" 
 res@cnLevels = (/ 0.1, 0.2, 0.5, 1, 5, 10, 15, 20, 25, 50/) 

 ;res@cnLevelSelectionMode =  "ManualLevels"
 ;res@cnMaxLevelValF       = 0.2
 ;res@cnLevelSpacingF      = 0.01 ;;02 ;;1

 nbin = 25

 PDF2 = pdfx(vxh(:,12,:,:),nbin,opt) 
 PDF3 = pdfx(vxh(:,13,:,:),nbin,opt) 
 PDF4 = pdfx(vxh(:,14,:,:),nbin,opt) 
 PDF5 = pdfx(vxh(:,15,:,:),nbin,opt) 
 PDF6 = pdfx(vxh(:,16,:,:),nbin,opt) 
 PDF7 = pdfx(vxh(:,17,:,:),nbin,opt) 
                    
 PDF = new((/8,nbin/),"double") 

 PDF!0 = "lev" 
 ;;PDF&lev = (/ 18, hgt(11), hgt(12), hgt(13), hgt(14), hgt(15), hgt(16), 7 /) 
 PDF&lev = (/ 18, hgt(12), hgt(13), hgt(14), hgt(15), hgt(16), hgt(17), 7 /) 

 PDF!1 = "vw" 
 PDF&vw = fspan(0.001,3.0,nbin) 

 print(PDF&lev) 

 PDF(0,:) = 0.
 PDF(1,:) = (/ PDF2 /) 
 PDF(2,:) = (/ PDF3 /) 
 PDF(3,:) = (/ PDF4 /) 
 PDF(4,:) = (/ PDF5 /) 
 PDF(5,:) = (/ PDF6 /) 
 PDF(6,:) = (/ PDF7 /) 
 PDF(7,:) = 0.

 plot2 = gsn_csm_contour(wks,PDF,res)

 overlay(plot(0),plot2) 

;***********************************************************
; panel resources
;***********************************************************
  panel_res                    = True   ; panel mods desired
  ;panel_res@gsnFrame           = False  ; don't advance frame yet
  ;panel_res@gsnPanelBottom     = 0.0    ; space for label bar
  ;panel_res@gsnPanelTop        = 0.5    ; only panel on lower half of page
  ;gsn_panel(wks, plot(0), (/ 1, 1 /), panel_res)

  gsn_panel(wks, plot(0), (/ 1, 1 /), panel_res)
  ;;gsn_panel(wks, plot2, (/ 1, 1 /), panel_res)




